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Power grids, road maps, and river streams are examples of infrastructural networks which are highly 
vulnerable to external perturbations. An abrupt local change of load (voltage, traffic density, or water level) 
might propagate in a cascading way and affect a significant fraction of the network. Almost discontinuous 
perturbations can be modeled by shock waves which can eventually interfere constructively and endanger 
the normal functionality of the infrastructure. We study their dynamics by solving the Burgers equation 
under random perturbations on several real and artificial directed graphs. Even for graphs with a narrow 
distribution of node properties (e.g., degree or betweenness), a steady state is reached exhibiting a 
heterogeneous load distribution, having a difference of one order of magnitude between the highest and 
average loads. Unexpectedly we find for the European power grid and for finite Watts- Strogatz networks a 
broad pronounced bimodal distribution for the loads. To identifj^ the most vulnerable nodes, we introduce 
the concept of node-basin size, a purely topological property which we show to be strongly correlated to the 
average load of a node. 



B 



lackouts, traffic gridlocks, and floods are all malfunctions of infrastructures which drastically affect their 
performance^ In many situations, they occur abruptly and might propagate through the network as shock 
waves^"^. These waves can either weaken by shedding their impact among branches or interfere construc- 
tively when two or more branches meet at the same node. The global consequences of these perturbations will 
strongly depend on the propagation dynamics and the capacity of each network element to bear abrupt 
changes^^"^^. The identification of vulnerable spots is a challenging scientific and technological question and this 
is precisely what we address here. 

Propagation of failures and cascading in complex networks have been subject of much scientific interest^^"^^ 
Examples are the use of the theory of self-organized criticality to study the propagation of failures in power grids 
and water transport on reservoir networks the Olami-Feder- Christ ensen model for earthquakes^^'^^, 
traffic^'^^'^^ and financial networks^^. Typically, the focus is on the cascading of failures resulting from an initial 
triggering event. However, it is also crucial to understand the dynamics preceeding these failures and identif)^ the 
vulnerable spots where they can possibly be triggered. 

To describe the propagation of shock waves on directed networks we use the Burgers equation^^. This equation 
describes flow when the flux depends quadratically on the load (e.g., voltage, traffic density, and water level). The 
range of applications of the Burgers equation goes beyond fluid dynamics as it is applied in many propagation 
processes, such as traffic jams, glacier avalanches or chemical processes^^'^^. Here we show that, in the case of 
perturbations randomly distributed in space, the dynamics of the solutions of the dissipative Burgers equation 
converges to a steady state in which the load distribution is strongly heterogeneous. Surprisingly, we find that the 
load of some nodes can exceed the average load by one order of magnitude. One might expect that the location of 
such nodes mainly depends on the propagation dynamics. Yet, we show that their fate is deeply imprinted in the 
network topology. We propose a new topological measure which allows to identif)^ the most vulnerable nodes 
without solving the dynamics. 



Model 

Dynamics. To describe the propagation of load (e.g., traffic density or water level) on a directed network, we 
consider on each link the one -dimensional Burgers equation^^ 



dp dp ^ 
dt ^ dx 



(1) 
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which we solve using Godunov's scheme. The details of the 
discretization and numerical solution are presented in the section 
Methods. 

Perturbation. Initially, the load on all directed edges and nodes is set 
to zero. Perturbations are described as local changes in the load 
according to the following procedure. First, we choose a node / at 
random and set p/ to a fixed value po (po > Pd during a time interval 
Tp. The load on the corresponding edges and on the other nodes is 
determined by solving the Burgers equation as described in the 
section Methods. After Tp, the constraint on the load of node / is 
released and its load is determined by the dynamics. A new node is 
selected and perturbed and the procedure is iterated. In addition to 
the perturbations, at each iteration step, 0.1% of every node load is 
dissipated. This dissipation would correspond, for example, to the 
evaporation of water from a river network, cars leaving the streets, or 
a potential drop due to Joule heating. 

Directed networks. The dynamics is investigated on the European 
high-voltage power grid^^ and two network models: the configu- 
ration model with power-law degree distribution^°"^^ and the 
Watts- St rogatz model, with small-world features^^. In the case of 
the model networks, the length of the edges are random variables 
chosen uniformly from the interval [3 : 20]. Initially, the power grid 
and the model networks are undirected. Inspired by the fact that in 
power grids the direction of the current depends on the node 
voltages, we use the following method to define the direction of the 
link. To each node /, a random value (the node voltage) is assigned 
uniformly from the interval [0:1] and the edge between two nodes is 
directed from the node with higher voltage to the one with lower 
voltage (i.e., (^source > target)- This method for generating directed 
edges automatically prevents the presence of loops. Since 
fluctuations are always present in the network, the direction of the 
current can vary in time. Our results are averaged over different 
voltage distributions as well. 

Results 

Steady state. At each time step, we measure the temporal load 
correlation, defined as: 



A-p{t) = 



p{t')-p{t'-mTp^ 



m 



(2) 



where is the load averaged over all nodes at time t. The brackets 
represent an average over the last ten consecutive time intervals of 
length lOOTp, i.e., t' = t - nlOOTp with n = 0, 1, 10. 

Starting with all loads equal to zero, we observe that Ap decays in 
time towards a steady state in which the dissipation balances the total 
incoming load. When Ap{t) drops below 1% of the relative standard 
deviation of the loads within the network, we assume that the steady 
state is reached. In the steady state, the load at each node has a well- 
defined average value with small fluctuations. The spatial distri- 
bution for a given realization of voltages in the European power grid 
is shown in Fig. 1. The size of the dots represents the average load 
over a time window of lOOOT^ measured in the steady state. Most 
nodes accumulate negligible load (p ~ 0), while surprisingly a small 
fraction of the nodes are overloaded (p > lOpo, where po is the 
magnitude of each perturbation). 

The observed load exemplarily shown in Fig. 1 corresponds to one 
realization of the voltage distribution and thus for one configuration 
of the direction of the edges. Assuming small temporal changes in the 
network (e.g., fluctuations of voltage in the power grid or number of 
cars entering a road junction), the direction of the edges changes in 
time. Thus, we also consider different realizations of the voltage 
distribution and, for each realization, we determine the steady state 
load distribution. Figure 2 shows the relative load distribution aver- 
aged over 5000 reaUzations. In order to compare the load distribution 




Figure 1 | Spatial distribution of load on the European power grid in the 
steady state. Distribution is shown for a particular realization of the 
voltage distribution. The size of the nodes corresponds to the average load 
allocated on them. The smallest size refers to zero load, where the largest 
one to the maximum load. 

of different networks (power grid, Watts-Strogatz and scale-free 
networks), loads in each curve are rescaled by the magnitude of each 
perturbation (po). The strongly inhomogeneous behavior of the 
steady-state load seen in Fig. 1 is also visible in the load distribution. 
The distributions are bimodal defining two different types of nodes: 
those with a negligible load compared to the perturbation (p < po) 
and those with a larger load (p > po). The latter ones are typically 
overloaded in the steady state, suggesting that the incoming pertur- 
bations interfere constructively at them. 

The plots for the two network models (Watts-Strogatz and scale- 
free) in Fig. 2 are obtained for networks with the same number of 
nodes as the power grid. The average degrees are also kept close to the 



o 
c 
0) 

CD 



10- 



2 10- 



CD 
> 

CD 
DC 



i5 10- 



10 



-12 





l| 1 1 1 1 1 IM| I| 1 






Power grid 




-m- Scale-free 




Small-world 





10- 



10-2 10-^ 10° 10^ 
Relative steady state load 



10^ 



Figure 2 | Relative load distribution in the steady state for different 
network topologies. The considered topologies are: power grid (green 
dots), scale-free (red squares) and Watts-Strogatz networks (blue 
triangles). The power grid and the scale-free network have N = 1254 nodes 
and M = 1811 edges, while the Watts-Strogatz network has N = 1254 
nodes and an average degree of (kout) = 2. All loads are in units of the 
amplitude of the perturbation. Each curve is an average over 5000 voltage 
realizations and, in the case of the model networks, also an average over 100 
different networks. The magnitude of the standard deviation of the curves 
is comparable to the size of the symbols. 
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power grid, with the same number of edges in the scale-free network 
and (kout) = 2m the Watts-Strogatz graph. In both cases, a bimodal 
distribution is also observed. The power-grid network is constructed 
from real data and its size corresponds to the real network size. Thus, 
a finite-size study is not possible. Yet, in the case of the model net- 
works one can systematically study the effect of the network size on 
the load distribution. Figure 3A shows the load distribution for 
Watts-Strogatz networks of different network sizes. The majority 
of the nodes (more than 90%) has always a negligible load, while 
the load of the remaining nodes follows a broad distribution, char- 
acterized by a decay in the relative frequency with increasing load 
and a cut-off for values of load close to po- The bimodal distribution 
smoothens out for larger network sizes. For scale-free networks the 
qualitative picture is slightly different. As shown in Fig. 3B, for all 
network sizes, one observes two power-law regimes, with a crossover 
at Pq. Nevertheless, note that for both network models, there is always 
a significant fraction of nodes (around 10%) with a non- negligible 
load. The load value of the cutoff suggests that, at large system sizes. 



A 




10-4 10-3 10-2 10-^ 10° 10^ 
Relative steady state load 



consecutive shock waves that enter the network are separated so that 
they attenuate their amplitude before being able to interfere. 

The specific nodes that exhibit these high load values typically 
change from realization to realization. However, after averaging over 
different voltage distributions, we still find some nodes which are 
consistently overloaded. For each distribution of voltages, we classify 
as "overloaded nodes" the ones with a load at least ten times larger 
than the average. We define vulnerability of a node as the probability 
that it is an overloaded node. Figure 4A shows the spatial distribution 
of vulnerability in the European power grid where the color and size 
of the nodes denotes their vulnerability. The vulnerability of white 
nodes is lower than 0.1%, while the one of the black nodes is larger 
than 5%. All the other nodes (about 30% of the nodes, in gray color) 
have a vulnerability between 0.1% and 5%. In comparison, the highly 
vulnerable nodes are at least 50 times more frequently exposed to 
large incoming fluxes. In the case of random perturbations, vulner- 
able nodes are more likely to fail or be congested. It is therefore 
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Figure 3 | Size-dependence of the relative load distribution in the steady 
state for two network models. The considered topologies are: (A) Watts- 
Strogatz and (B) scale-free networks with (k) = 2. Insets show the 
respective data collapse, where p and P denote the load and the relative 
frequency, respectively, N is the size of the network. Loads are divided by 
the magnitude of the applied perturbation to ensure comparability. Each 
data is averaged over at least 100 different graphs and 100 voltage 
realizations. The breaks in the distributions around p = 10"^ are due to the 
logarithmic binning. 




Figure 4 | Spatial distribution of two properties on the power grid 
network. (A) Vulnerability of the nodes obtained solving the Burgers 
equation (probability of having lOp in a voltage realization). (B) Node- 
basin sizes averaged over different voltage realizations. Color and size 
indicate the strength of the corresponding property: black corresponds to 
large vulnerabilities (> 5%) and basin sizes (> 10), while white nodes have 
negligible vulnerabilities (below 0.5%) and small basin sizes (close to 1). 
The data are the average over 5000 voltage realizations. 
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Figure 5 | Distribution of vulnerabilities for three different network 
topologies. The distribution of the vulnerability (i.e., the probability of 
having a load ten times larger than the average) is shown for power grid 
(green circles), scale-free network (red squares) and Watts-Strogatz graph 
(blue triangles). The data for the model networks are the average over 100 
network realizations. 

crucial to identify these nodes to improve their capacity and to mit- 
igate the risk of failure. Figure 5 shows the vulnerability distribution 
corresponding to the map in Fig. 4A and to other network topologies. 
In the case of the Watts-Strogatz and scale-free networks, a vulner- 
ability distribution for networks of size AT = 10^ are presented. 

Identifying vulnerable nodes. Next we will introduce a simple 
topological property of the nodes to identify the vulnerable spots 
without solving the dynamics of the Burgers equation. According 
to the Burgers equation, each node sheds the incoming shock wave 
among its out-edges and the total load agregated at the node at a time 
t is the sum of incoming loads. Hence, if we track the path of a given 
shock wave, it is fragmented at each node with multiple out-edges 
and will stop at any node that does not have any out-edges. Assuming 
that the time average of the load at a node is proportional to the 
number of incoming shock waves, we determine the basin 
corresponding to that node. For this purpose, let us consider one 
realization of the edge directions (see Fig. 6). For a given node (the 
one in red in Fig. 6), the corresponding basin is defined as 
the smallest subgraph of the network containing this node (as the 
sink), which is connected to the rest of the network by outgoing 
edges. Following the procedure as illustrated in Fig. 6, we go 
through the nodes following the opposite direction of the in- edges, 
starting from the red node (Fig. 6A), and add nodes to the basin that 
has only out-edges at the end of the process (marked by the blue 
region in Fig. 6B). The resulting subgraph is the basin of the red node. 



and the contribution of the load of the nodes in the basin is simply the 
inverse of their out-degree (Fig. 6C), except for the initial (red) node, 
that contributes to the load with unity. This choice of the 
contribution of the nodes in the basin is based on the fact that Eq. 
(7) conserves the flux at each node. For simplicity, we assume here 
that the amplitude of the Shockwaves leaving a node is on average the 
same for each out- edge. 

We calculated the size of the basins for each node, defined as the 
sum of contributions from all nodes inside the basin of the corres- 
ponding node. This basin size (which is determined for one voltage 
distribution) is then averaged over different voltage configurations. 
The resulting basin size distribution is depicted in Fig. 4B for the 
power-grid network. One sees that for this network, the distribution 
of the node-basin size is very similar to the distribution of the vul- 
nerability. A quantitative comparison of the two properties can be 
given by their correlation. Thus, we plot the rank- rank scatter plots in 
Fig. 7 of the indices of the nodes after being sorted in ascending order 
by vulnerability and node-basin size. The plots are the average over 
5000 voltage distributions and over 100 different topological realiza- 
tions of model networks. 

The corresponding product-moment correlations p^^ of the ranks 
are given above the plots, showing strong correspondence between 
the ranked vulnerability and node-basin size for the power grid and 
scale-free networks. The product-moment correlation is. 



(3) 



where fi^ and fiy are the mean and and Gy are the corresponding 
standard deviation of the two quantities x and y. The crucial nodes 
are the ones with large vulnerability since they are more exposed to 
large loads. The basin size shows a strong correlation with the vul- 
nerability for these large values (top right corner of the scatter plots), 
meaning that it is a good estimator of vulnerability. Watts-Strogatz 
network exhibits much less correlation because the degree distri- 
bution is extremely narrow, that is, deviations from the average 
degree are negligible. Thus, for each realization, the differences in 
the loads from node to node are very small. This conclusion is also 
supported by Fig. 5 showing a narrow vulnerability distribution for 
Watts-Strogatz networks. In order to describe the nodes that are 
overloaded more frequently than the average in the small-world 
network, we measured the average length of the in-edges of the 
nodes, i.e., the average number of mesh points on the in-edges. 
Also, to obtain a more realistic distribution of the edge lengths, the 
graphs were constructed by first considering a two-dimensional lat- 
tice and then rewiring each edge with probability p = 0.01, finally we 
applied a deviation of 50% to each length, distributed uniformly 
around the original values. Simulations on 50 network realizations 
with AT = 32 X 32 nodes and 1000 voltage realizations for each 
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Figure 6 | Calculation of the node-basin size. (A) The basin corresponding to the red node is considered. (B) We determine the smallest subgraph 
containing the red node and having only out-edges to the rest of the network. This is equivalent to a breadth-first search traversing the fraction of the 
graph reached through only in-edges. (C) When the basin is determined, each node in the basin contributes to the red node's basin size by the inverse of its 
out-degree. The contribution of the red node (i.e., the sink of its basin) is unity. 
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A p^^ = 0.95 B p^^ = 0.93 ± 0.006 C p^^ = 0.36 ± 0.034 




Vulnerability Vulnerability Vulnerability 



Figure 7 | Scatter plots of the ranked vulnerability and node-basin sizes for different network topologies. The corresponding network models are: (A) 
power grid, (B) scale-free and (C) Watts-Strogatz network. The axes on the plots indicate the corresponding ranked property, and the color denotes the 
number of nodes for which the two ranks were identical. In other words, the color of a dot at (x, y) corresponds to the number of nodes with 
vulnerability rank of x and node-basin size rank of y (see the colorbars). All plots are the average over 5000 voltage realizations. In the case of model 
networks, 100 different realizations are considered. 



network resulted in a Pearson-correlation of -0.51 between the vul- 
nerability of the nodes and the average length of their in-edges (with 
Pearson-correlation of 0.1 between vulnerability and basin size), 
meaning that nodes with short in-edges are more frequently over- 
loaded. On short edges, decrease in the magnitude of the Shockwaves 
is lower, compared to the attenuation of the waves on longer ones. 
This suggests that in the case of small-world networks, structural 
homogeneity is very strong (differences among nodes are less pro- 
nounced, compared to the scale-free networks), and the physical 
length of the edges through which load is transmitted have to be 
considered. See Supplementary Material for snapshots and the scat- 
ter plot of the average length and vulnerability in this small-world 
network model. 

Discussion 

We study the propagation of Shockwaves on directed networks using 
the Burgers equation. Under sequentially applied perturbations and 
constant dissipation, the dynamics approaches a steady state. In this 
steady state, most of the nodes have negligible average load and a 
significant fraction of the total load is localized on a few nodes. We 
found that some nodes are more likely to accumulate load even after 
averaging over many edge direction configurations. These nodes (the 
vulnerable nodes) are more likely to fail, when there is a finite capa- 
city of the load they can bear. Unexpectedly we find for the European 
power grid a broad pronounced bimodal distribution for the loads, 
while for scale-free network the distribution resembles more a power 
law. 

The steady state and thus the probability distribution of vulner- 
ability among the network is determined by solving numerically the 
partial differential equations of Eq. (1) on each edge. The propaga- 
tion velocity of the shock waves depends on their amplitude, which 
can vary rapidly throughout the network. We propose a simpler 
alternative based on the node-basin size to estimate the vulnerability 
of the nodes and identify the most vulnerable ones. Simulations on a 
real network (European high-voltage power grid) and on scale-free 
networks show that the node-basin size can predict very accurately 
the location of vulnerable nodes while it performs worse for the 
Watt-Strogatz network due to its narrow degree distribution. 
Investigation of the edge lengths suggests that in the Watts- 
Strogatz network, vulnerable nodes are related to those having very 
short in-edges. 

Our results suggest that it is possible to establish a remarkable 
connection between dynamics and network structure. Although for 
many networks the node-basin size seems to be an accurate tool in 
predicting the distributions and detecting vulnerability, it is only the 
first step towards a complete description of the steady state. More 
might be understood by studying the relation between the most 
vulnerable nodes: under what circumstances are they separated or 



forming connected subgraphs? Is any local property of the network 
responsible for a node being highly vulnerable? This information 
would provide the tools to mitigate the risk of systemic failure. 
Further investigation may involve the removal of nodes that reach 
their capacity. In this case, the study of the time evolution of the 
network structure or optimal strategies of dynamical node/edge 
addition or deletion can be of relevance. 

Methods 

Dynamics. In this section, we describe the generalization of Eq. (1) on a directed 
network. The numerical solution of the one- dimensional Burgers equation can be 
discretized using Godunov's scheme^^'^^ 

in — flow out — flow 

pr=p'i+^[HZ^-HM^], (4) 

where p- is the load at the mesh point / at time t. Ax and At are the spatial and 

temporal discretizations, and F{p,ri) = — is the flux. The value of p is given as 
follows^^: If p > ?7 then 

(3) 

I rj otherwise 

otherwise, 

IP if p>0 
rj if rj<0 . (6) 
0 if p<0<r] 

To solve this equation on a network, one needs to fix the direction of each link in order 
to have a precise definition of the in- and out-flux of a mesh point. For practical 
purposes, this is a realistic approach as water always flows downhill and the current 
follows a decreasing gradient in electric potential. Our discretization model for the 
edges and nodes is illustrated in Fig. 8. The edges of a network are one dimensional 
and thereby Eq. (4) holds. The number of mesh points in each edge is proportional to 
its length and the direction is defined by the direction of the edge. In the case of model 
networks, the length of the links is a uniformly distributed random value. 
Furthermore, a value of pk is assigned to each node k. Nodes interact with the nearest 
mesh points of their incident edges according to the following equation, 

in — edges 
out — edges 

gin f^gout^ denotes the set of in- (out-) edges of node /, and (r^) is the load at the last 
(first) mesh point of the corresponding edge. The resulting dynamics conserves the 
total mass and at each node, the total incoming flux is equal to the outgoing one. 
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order to eliminate degree-correlations. Note that the number of edges in the Watts- 
Strogatz network is different from that in the power grid and the scale-free network. 
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Figure 8 | Illustration of the discretization model on a directed network. 

(A) Each edge corresponds to a one-dimensional coordinate system with 
the positive direction defined by the direction of the edge. Xi denotes the rth 
mesh point on the edge, Xs and x^are the mesh points corresponding to the 
source and target nodes, respectively. The load values on the ith edge mesh 
point are denoted by fj. (B) The nodes interact only with the nearest mesh 
points of their edges, r/ denotes the adjacent mesh point of the 
corresponding in- or out-edge of the blue node and Fis the flux according 
to Eq. (4). 



Also, we should add an important remark on the various constraints of the 
investigated model. First, in the numerical solution of a nonlinear PDE on a network, 
the degree of a node corresponds to the local dimension of the space in which one 
solves the equations. As the size of the scale-free network increases, the frequency of 
nodes with very large degrees also grow. Considering numerical stability, the 
appearance of larger degrees sets an upper limit on the magnitude of the applied 
perturbations. However, if the perturbations are small (which is required by the 
numerical treatment), shock waves tend to vanish by travelling on the edges and they 
are not able to interfere constructively. Therefore, in the finite- size study, we con- 
sidered only networks below the size of AT = 10^. 

Networks. We consider three network models: the European power grid (with N = 
1254 nodes and M = 1811 edges, i.e., (kout) = 1-44), the Watts-Strogatz small-world 
network {{kout) = 2) and scale-free network models {(kout) = 1-44 for N = 1254 and 
(kout) = 2 in the case of large network sizes). The Watts-Strogatz network is 
constructed by considering first a one- dimensional chain with first and second 
neighborhood connections and periodic boundary conditions, and then rewiring 
each edge with probability p = 0.01 (with undirected edges, this corresponds to an 
undirected average degree {k) = 4). After the voltages are set and edge directions are 
introduced, the resulting network has a directed average degree of (kout) = 2. The 
scale-free network is constructed by the configuration model: first we assign the 
degrees for each node according to a power-law with exponent y = 2.5, and then 
connect randomly chosen nodes. Finally, further rewiring of the edges is carried out in 
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